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THE  ITERATED  SIMPSON  METHOD 
CF  NUMERICAL  INTEGRATION 

ABSTRACT 

The  paper  reviews  the  well  known  Iterated  Simpson  method  for 
obtaining  the  numerical  solution  of  the  Initial  value  problem  for 
a  general  system  of  n  first  order  ordinary  differential  equations® 
The  formulae  are  developed  and  two  examples  are  given  illustrating 
their  use® 
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I,  INTRODUCTION 


Earl7  in  1950  a  new  numerical  procedure  for  solving  systens  of 
ordinary  differential  equations  on  conputing  machines  was  announced 
by  Rs  Ft,  Clippinger  and  Bo  Dimsdale  of  the  Computing  Laboratory,,  The 
method  was  first  discussed  in  the  Coding  Notes  prepared  for  the  Coding 
Course  given  by  the  Ballistic  Institute  (September  19h9  to  February 
1950)®  The  procedure  is  also  mentioned  as  Clippinger ’'s  Method  in  the 
lecture  notes  of  L®  Ho  Thonas®  Ifowever^,  since  neither  set  of  notes  is 
available  for  general  distribution  it  was  felt  desirable  to  publish 
this  note. 

This  method^  which  has  been  called  "Iterated  Siupson"^  is  of  the 
foxirth  orders  that  is  the  error  (excluding  round-off)  goes  to  zero  as 
the  fifth  power  of  the  step  size®  It  is  thus  equivalent  to  approxi¬ 
mating  the  solutions  locally  (on  intervals  of  leixgth  2h)  by  polynomials 
of  the  fourth  degree®  These  polynomials  differ  from  the  Taylor  expan¬ 
sions  of  the  dependent  variables  about  the  left  endpoints  of  the  inter¬ 
val  of  integration  by  terms  of  order  only 5  if  the  right  members  of 

the  differential  equations  possess  a  sufficient  number  of  bounded  deri¬ 
vatives  » 

like  the  Runge-Kutta  method^  it  proceeds  from  initial  data  at  one 
point  without  need  for  data  at  other  points®  This  is  a  distinct  advan¬ 
tage  for  large  scale  digital  computers  since  it  limits  the  storage  re¬ 
quired  for  the  dependent  variables ^  and  simplifies  changes  in  the  grid 
size  idien  required®  If  the  functions  possess  discontinuities  j,  they 
may  be  placed  at  the  grid  points  and  the  method  may  still  be  used® 

Another  advantage  is  the  smaller  amount  of  conputer  memory  required 
to  solve  a  system  of  given  order,  when  conpared  with  other  methods  of 
conparable  accuracy®  It  has  been  claimed  that  systems  of  roughly  60% 
higher  order  can  be  solved  with  this  method  than  with  the  widely  used 
methods  of  Moulton,  HLlne,  Adams,  etc®  However,  the  method  of  Runge- 
Kutta  can  be  shown  to  be  as  efficient  as  Iterated  Sinpson  with  respect 
to  storage  requirements,  for  Instance,  Proceedings  of  the  Cambridge 
Philosophical  Society,  Vol®  h7  Biirt  1  1950,  "A  Process  for  the  Step  by 
Step  Integration  of  Differential  Equations  in  an  Automatic  Digital  CoI1^'^ 
putlng  Machine"®  S®  Qill,  pp  96-109®  Further  atud^  will  indicate 
idiether  the  method  of  Runge-Kutta  as  used  by  S®  Qill  should  replace 
Iterated  Sinpson  entirely® 

Further  it  should  be  enphasieed  that  the  usefulness  of  the  Iterated 
Sinpson  method  is  restricted  largely  to  high  speed  digital  conputers, 
since  considerable  conputatlon  time  is  sacrificed  in  order  to  gain 
these  advantages®  However,  for  a  fast  machine  with  limited  storage 
this  saorlfloe  may  well  be  oonpensated  for  ^  the  smaller  storage  re¬ 
quired®  A  distljiot  advantage  over  Runge-Kutta  is  that  the  Iteration 
makes  a  single  step  practically  self  oheoking® 
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IIo  PROCEDURE 


Let  the  system  of  n  first  order  ordinary  differential  equations  be 

(1) 

idiich  can  be  written  in  vector  form  as  y^  ■  f (x^iy)  idaere  y,y’  and  f 
represent  vectors.  The  initial  conditions  are 

(2)  yCx^)  “  yi„i  at  X  - 

He  derive  fox^mulae  for  y  at  a  point  x  which  is  a  distance  2h  from 
by  expanding  y(x)  »  y^s  ya(x)  »  yj^j  y(x^)  »  and  y'(x^j)  -  y^^^ 

about  X  ®  ^  and  eliminating  y  ®  yCx),  y%  y",  y"*  j,  and  y*'''"  from 

the  expans ions s 

{3«1)  y^^  “  y  ■>■  hy«  ♦  (h^/2)  y«  4-  (h^/6)  y“®  +  (h^Sl;)  y*"^  +  O(h^) 

(3.2)  y^^^^  ffl  y  -  hy«  (h^/2)  y«  =  (l?/S)  +  (hVsU)  y*"^  *  (h^) 

(3.3)  y|  -  y»  +  h^  +  ih^/2)  +  (h^/e)  y"^  O(h^) 

(3.1;)  y“^^  »  yo  -  hy"  +  (hVz)  y"»  »  (h^/6)  y'V  +  o(h^), 

where  0(h™)/h’’^  is  a  function  of  h  which  renains  bounded  as  h  -*•  0, 
Adding  (3,1)  and  (3.2)  we  get 
(U.l)  y^^  +  =  2y  +  h^^i  +  O(h^) 

Subtracting  (3,2)  from  (3.1)  ve  get 
(U.2)  y^^  »  yj^_^  “  2hy»  +  h^/3  y""  +  O(h^) 

Adding  (3.3)  and  (3.U)  we  get 

(1;.3)  y'  y®  »  2y«  +  h^p*'  +  O(h^) 

^  4-“l 

Subtracting  (3.U)  from  (3.3)  we  get 

(U.U)  y«  “  y !  _  “  Zhy"  ♦  hV3  y'"^  *  o(h^), 

i  i“l 

Next  we  eliminate  y"®  by  multiplying  (l;.3)  by  -h/3  and  adding  it  to  (1;.2) 
Collecting  terms  we  get 
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(5.1)  7^  ®  7^^]^  +  ^  4-  1;7«  +  7^)  O(h^) 

Equation  (5.1)  is  iiierel7  Sijjp3on''s  Rule  which  is  knovm  to  be  correct 
to  fourth  order® 

Likewise  we  eliminate  7”  b7  multiplTing  (U.l)  b7  2  and  adding  to 
“  h  times  (U.U).  Divide  this  relation  through  b7  h  and  we  get 

(5.2)  7  ~  7^  *  7^^i  ^^i-1  ”  ^i^  ^  0(h**) 

The  numerical  procedure  for  determining  7(1)  consists  in  taking^, 
as  a  first  approximation  b7  Euler^ 


(6) 
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"i=l 


and  defining  the  (j  +  l)st  approximation  to  7^  b7  the  formulaes 


(7.1)  7j  «  7ij  /2  +  7i„i  /2  +  hA  (7[^^  =  7|^j)  and 

■  ^i-i  ♦  *>■'•3  (yu  ♦  *^3  ♦  yij>” 

^1  7®  in  |7«1)  and  (7.2)  are  defined  hy  7'  »  f(xg7)s  in  particular 
7'  ®  ^(5.9  7j). 

Before  listing  a  few  exanples  we  euimnarlze  how  the  method  is  needs 

a.  First  use  Euler  integration  to  get  an  initial  approximation  to 
the  values  of  the  dependent  variables  at  the  new  point. 

bo  Next  use  these  values  and  the  Talues  at  the  initial  point  to 
obtain  approximations  for  the  values  of  the  dependent  variables 
at  the  point  half  uay  between  the  new  point  and  the  old  point. 

Ca  Use  these  intermediate  values  and  the  values  at  the  old  point 
to  obtain  a  better  approximation  to  the  new  point. 

d.  Steps  b.  and  c.  ]na7  be  repeated  as  inan7  times  as  desired  or 
until  the  changes  in  the  dependent  varLablea  are  sufficientl7 
small. 


In  the  exan^lee  considered  in  the  next  section two  Iterations 
were  found  to  be  sufficient  to  obtain  the  desired  accuracy.  In  other 
exan^jles  more  than  two  iterations  were  required.  It  should  be  noted 
that  the  ordering  of  the  steps  is  such  that  only  the  last  confuted 
values  at  the  new  point  or  the  values  at  the  intemsdiate  point  need 
be  kept  in  storage  in  addition  to  the  values  at  the  old  point. 
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m.  EXAMPLES 


Us  illtistrate  the  procedure  b7  means  of  tiio  examples |  both  were 
coded  for  and  run  on  the  Eniac. 


as  Bessel  Functions  of  Order  zero  and  one. 


The  method  ^s  first  applied  to  the  solution  of  the  equations 

(8*1)  y*  “  “  z  and  (8,2)  z‘  «  y  -  z/x  with  initial  conditions? 

(8,3)  y(0)  ■■  Ip  z{0)  0,  ;The  solution  of  this  system  is 

y  «  Jq(3c)  (zero  order  Bessel  function) 
a  “  (first  order  Bessel  function) 

It  was  found  possible  to  obtain  nine  correct  significant 
figures,  comparing  the  results  with  the  Harvard  Tables,  The 
error  due  to  round-off  was  only  a  maximom  of  5  in  the  10th  figure 
after  500  steps  of  integration  (an  estimated  30,000  multiplications 
and  divisions).  This  is  only  10  times  the  error  involved  in  intro¬ 
ducing  a  number  into  the  Eniac, 

The  algorithm  for  Iterated  Slnpson  applied  to  the  problem 
(8,1)  -  (8.3)  and  used  on  the  Eniac  follows? 

As  before  we  used  (6)  for  the  initial  approximations  (or  for 
the  two  equations  we  write). 


(8,U)  y^^^,  =  y^^^  +  ^^i-1 

(8,5)  *  ^*“l-l 

The  Iterated  Simpson  relations  are 


(8.6)  2/3  y^j  "  y^j  /3  +  7^^^  /3  +  h/6  (z^j  - 

(8.7)  2/3  z^j  a  z^^  /3  +  /3  +  h/6  “  (f^j  "  /x)^ 

<h-l  *  ‘^±3  * 

*1,3,1  ■  *1-1  ♦  +  1*  <yij  -  /S)  ♦  yi,  -  A 

2/3  y^j  is  computed  rather  than  y^^  to  aid  in  handling  number  size 
and  to  sliiplify  the  number  of  operations  required. 


b,  Integration  of  the  trajectory  equations. 
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The  method  has  been  applied  to  the  integration  of  the  particle 
theory  trajectozy  eqmtions  uith  t  as  the  independent  variable. 

The  algorithm  follows: 

For  the  initial  approxination  we  warite 

(9.1)  y^^  -  y^_^  “  ^®i-l^i„l  ♦  g)  2h 

(9.2)  Zh 

(9.3)  =  Ti.i  ♦  21iy^ 

(9.U) 


The  Iterated  Simpson  relations  are  as  follows: 


(9.5) 

*  * 

hA 

H 

9 

(9.6) 

X®  “■ 

ij 

1/2  (x.® 

i”! 

*  " 

hA 

”  =i- 

(9.7) 

y.  .  " 

1/2  (y^^]^ 

*  5^1^)  * 

hA 

- 

y*  ) 

^±y 

(9.8) 

a 

^ij  “ 

l/2(x^_^ 

*  ♦ 

hA 

- 

(9.9) 

yO 

»  y  ®  = 

^i-1 

'  h/3  (E., 

-in. 

■1  *  ‘^13^13  * 

=13=^13 

(9.10) 

.  "  7i-l  * 

h/3 

♦ 

■1 

(9-11) 

x' 

“  ^-1  “ 

h/3  (E^, 

“l^i- 

.1  *  ‘^13^i3  * 

*13*13^ 

(9.12) 

“  *i-l  ♦ 

h/3  (x® 
1- 

■1  * 

ijx®^  ♦  3 
ij 

c®  ) 

A  discussion  on  bounds  on  errors  is  contained  in  the  original 
Coding  Notes  of  CUpplnger  and  Cimscjale^  as  well  as  in  the  Lecture 
Notes  on  "Ordinaiy  Differential  Bquationsg  June  -  August  1950" »  pre=> 
pared  by  Iferk  Lotkin^  which  are  available  oiily  for  local  BEL  (3istri“ 
but ion. 


Ttf.  Barkley  Fritz 
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